/************************************************************************
 * MechSys - Open Library for Mechanical Systems                        *
 * Copyright (C) 2013 Sergio Torres                                     *
 *                                                                      *
 * This program is free software: you can redistribute it and/or modify *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * any later version.                                                   *
 *                                                                      *
 * This program is distributed in the hope that it will be useful,      *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the         *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with this program. If not, see <http://www.gnu.org/licenses/>  *
 ************************************************************************/

// STL
#include <iostream>
#include <cmath>

// MechSys
#include <mechsys/util/fatal.h>
#include <mechsys/util/array.h>
#include <mechsys/util/util.h>
#include <mechsys/linalg/matvec.h>
#include <mechsys/mesh/structured.h>
#include <mechsys/mesh/unstructured.h>
#include <mechsys/dem/domain.h>

using std::cout;
using std::endl;
using std::ifstream;
using std::min;
using std::max;

void GenFromMesh(DEM::Domain & dom, char const * VertsFile, char const * CellsFile, double thickness, double scale, double R, double rho)
{
    String fnv(VertsFile); fnv.append(".mesh");
    String fnf(CellsFile); fnf.append(".mesh");
    ifstream fnvf(fnv.CStr());
    ifstream fnff(fnf.CStr());
    if (!Util::FileExists(fnv)) throw new Fatal("File <%s> not found",fnv.CStr());
    if (!Util::FileExists(fnf)) throw new Fatal("File <%s> not found",fnf.CStr());
    size_t ncol = 0;
    Vec3_t Vtemp = OrthoSys::O;
    size_t number;
    Array<Vec3_t> Vtotal;
    while (!fnvf.eof())
    {
        if (ncol==0) fnvf >> number;
        if (ncol==1) fnvf >> Vtemp(ncol-1);
        if (ncol==2) fnvf >> Vtemp(ncol-1);
        ncol++;
        ncol = ncol%3;
        if (ncol==0)
        {
            Vtotal.Push(scale*Vtemp);
            //std::cout << Vtotal[Vtotal.Size()-1] << std::endl;
        }
    }
    Array<Array <size_t> > Ctotal;
    Array<size_t>        Cpartial;
    size_t ncoltot = 0;
    ncol = 0;
    while (!fnff.eof())
    {
        if (ncol==0)
        {
            fnff >> ncoltot;
            ncoltot++;
            Cpartial.Resize(0);
            Cpartial.Resize(ncoltot-1);
        }
        else
        {
            fnff >> number;
            Cpartial[ncol-1] = number;
        }
        ncol++;
        ncol = ncol%ncoltot;
        if (ncol==0)
        {
            Ctotal.Push(Cpartial);
            //std::cout << Ctotal[Ctotal.Size()-1][0] << " " << Ctotal[Ctotal.Size()-1][1] << " " << Ctotal[Ctotal.Size()-1][2] << " " << Ctotal[Ctotal.Size()-1][3] << " " << std::endl;
        }
    }

    //for (size_t i=0;i<Ctotal.Size();i++)
    //{
        //for (size_t j=0;j<Ctotal[i].Size();j++)
        //{
            //std::cout << Ctotal[i][j] << " ";
        //}
        //std::cout << std::endl;
    //}

    //return;


    Array<Vec3_t> V;
    Array<Array <int> > E;
    Array<Array <int> > F;
    for (size_t i=0;i<Ctotal.Size();i++)
    {
        size_t nverts = Ctotal[i].Size();
        V.Resize(2*nverts);
        E.Resize(3*nverts);
        F.Resize(nverts+2);
        //std::cout << "1" << std::endl;
        for (size_t j=0;j<Ctotal[i].Size();j++)
        {
            V[j+nverts] = Vtotal[Ctotal[i][j]] + Vec3_t(0.0,0.0,-thickness);
            V[j]        = Vtotal[Ctotal[i][j]];
        }
        //std::cout << "2" << std::endl;
        for (size_t j=0;j<nverts;j++)
        {
            E[j].Push(j);
            E[j].Push((j+1)%nverts);
            E[j+nverts].Push(j);
            E[j+nverts].Push(j+nverts);
            E[j+2*nverts].Push(j+nverts);
            E[j+2*nverts].Push((j+1)%nverts+nverts);
        }
        //std::cout << "3" << std::endl;
        for (size_t j=0;j<nverts;j++)
        {
            F[j].Push(j);
            F[j].Push(j+nverts);
            F[j].Push((j+1)%nverts+nverts);
            F[j].Push((j+1)%nverts);
            F[nverts].Push(j);
            F[nverts+1].Push(2*nverts-j-1);
        }
        //std::cout << i << std::endl;
        //for (size_t i=0;i<V.Size();i++)
        //{
            //std::cout << V[i] << std::endl;
        //}
        //for (size_t i=0;i<E.Size();i++)
        //{
            //std::cout << E[i] << std::endl;
        //}
        //for (size_t i=0;i<F.Size();i++)
        //{
            //std::cout << F[i] << std::endl;
        //}
        DEM::Erosion(V,E,F,R);
        dom.Particles.Push (new DEM::Particle(-1, V,E,F,OrthoSys::O,OrthoSys::O,R,rho));
    }
}

int main(int argc, char **argv) try
{
    String cellsfile (argv[1]);
    String vertsfile (argv[2]);
    DEM::Domain dom;
    GenFromMesh(dom,vertsfile.CStr(),cellsfile.CStr(),1.0,10.0,0.1,1.0);
    dom.WriteXDMF("dem_frack");
}
MECHSYS_CATCH


